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Abstract 

We analyze a two dimensional viscoresistive magneto- 
hydrodynamical (MHD) model for a thin accretion disk 
which reconciles the crystalline structure outlined in 
Coppi(2005)[ Coppi and Rousseau(2006)"] with real micro- 



pact object. Starting from the original idea of Shakura 
Shakura(1973)| a wide number of results (see for istance 



scopic and macroscopic features of astrophysical accreting 
systems. In particular, we consider small dissipative ef- 
fects (viscosity and resistivity, characterized by a magnetic 
Prandtl number of order unity), poloidal matter fluxes and 
a toroidal component of the magnetic field. These new 
ingredients allow us to set up the full equilibrium profile 
including the azimuthal component of the momentum con- 
servation equation and the electron force balance relation. 
These two additional equations, which were identically sat- 
isfied in the original model, permit us to deal with non 
zero radial and vertical matter fluxes, and the solution we 
construct for the global equilibrium system provides a full 
description of the radial and vertical dependence within 
the plasma disk. The main issue of our analysis is out- 
lining a modulation of the matter distribution in the disk 
which corresponds to the formation of a ring-like sequence, 
here associated with a corresponding radial oscillation of 
the matter flux. 



1 Introduction 

The relevance of the accretion processes in understanding 
the emission properties of many astrophysical sources like 
Gamma Ray Bursts [Piran(1999)| and Active Galactic 
Nuclei |Krolik(1999)1 , led to a significant effort over 
the last four decades in order to construct a standard 
model for the profile of an accreting disk over a com- 



Pringle and Rees(1972)|, [Shakura and Sunyaev(1988) , 



or |Lynden-Ben and Pringle(1974)[ |Lynden-Bell(1996)| ) 
allowed to achieve the firm understanding that, 
within a disk, the angular momentum transport is 
driven by dissipative effects (for a detailed review see 
Bisnovatyi-Kogan and Lovelace(200T)] ). Since the nature 



of such viscosity and resistivity of the plasma cannot 
be recovered on a microscopical level (such systems 
are mainly collinsionless) , then the origin of non-ideal 
features must rely on the turbulent behaviour of the 
plasma once linear instabilities are triggered. The most 
successful paradigm to explain the onset of the turbulence 
is the so called Magneto-Rotational Instability, discov- 
ered by Velikhov [Velikhov(1959)] and Chandrasekhar 
[Chandrasekhar(1960)| , and recently fully developed by 
Balbus and collab orators [Balbus and Hawley(1991)| (for 
a review see also [Balbus and Hawley(1998)| ~ However, 
such instability, which can be triggered by an arbitrary 
small magnetic field, is nonetheless suppressed when 
the sound speed is smaller than the Alfven velocity, i.e. 
when the /? parameter of the plasma is smaller than one. 
Since this condition can take place in sufficiently cold 
and magnetized plasma disks, a challenging puzzle arises 
when dealing with the accretion mechanism in structures 
exhibiting such regimes. 

Two interesting proposals for alternative scenarios 
were proposed in [Bisnovatyi-Kogan and Lovelace(2000)[ 



and Coppi and Coppi(2001)[ . In particular, B. Coppi 
and coworkers developed over the last ten years a line 
of research devoted to construct a different paradigm 
Coppi(2009)| , mainly based on the implementation on 



the astrophysical setting of some basic features observed 
in laboratory plasma, like ballooning mode instabilities 
Coppi and Keyes(2003)| . More than solving the problem 



of accretion for an ideal plasma disk, such efforts had the 
main merit to outline the existence of regimes where the 
disk profile can fragment into a sequence of rings. This sit- 
uation takes place as soon as the Lorentz force, produced 
by the plasma backreaction in the magnetic field of the 
central object, plays a significant role in the vertical con- 
finement. This result, however, was obtained in a local 
description of a purely rotating disk, and disregarding any 
dissipative effect (see also [Montani and Carlevaro(2010)| 
for the jet formation in this scheme). Despite being rather 
small, such dissipation features are present in any real as- 
trophysical system, and the possibility to reconcile the ring 
sequence with a realistic equilibrium configuration is the 
main task addressed by the present work. 

We consider the full MHD equilibrium including, with re- 
spect to the original model, the azimuthal momentum con- 
servation equation and the electron force balance relation. 
For a purely ideal and rotating plasma, these equations are 
identically satisfied, but in a real system (here we consider 
non zero viscosity and resistivity coefficient, associated to 
a magnetic Prandtl number of order unity) they are cru- 
cial for fixing the poloidal matter fluxes and the toroidal 
magnetic field. We provide a solution for the global con- 
figuration system which has the merit to demonstrate that 
the existence of the ring sequence within the disk is com- 
patible with non zero dissipative effects, and, over all, with 
the presence of a non zero radial matter fiux. Nonethe- 
less, the smallness of the additional terms we include in 
our treatment, allows us to neglect the poloidal velocity 
and the toroidal magnetic field in the radial and vertical 
configuration, which therefore retain the same form as in 
Coppi(2005)[ Coppi and Rousseau(2006)] . We emphasize 



that dealing with only small corrections is a constraint im- 
posed by the corotation theorem [Ferraro(1937)| . In fact, 
as shown in the last section, up to the leading order, it is 
still possible to preserve the relation between the angular 
velocity of the disk u and the magnetic surface function ?/; 
as far as the additional contributions due to the poloidal 
velocity, the toroidal magnetic field and the dissipative ef- 
fects, are sufficiently small with respect to the basic model 
quantities as in [Coppi(2005) Coppi and Rousseau(2006)"] . 
A significant improvement that our analysis is able to pro- 
vide, is fixing the vertical dependence of the disk profile 
as a consequence of a compatibility equation coming out 
from the azimuthal and the electron force balance system. 
The local model we provide can be regarded as the starting 
point for the study of the settlement of linear instabilities 
in the paradigm of the crystalline structure. 



As outlined in Montani and Benini(201l")] , the crys- 



talline morphology of the magnetic field, as well as the ring 
sequence profile, have a characteristic scale which turns out 
to be microscopical on an astrophysical setting (indeed the 



typical length of the perturbations results to be of the order 
of some meters). This issue suggests that, more than for 
its stationary appearance, such a morphology concerns the 
stability structure of the disk, and only after a full scheme 
of perturbation analysis is developed, a precise nature of 
the configuration can be settled down. In this respect, our 
treatment stands as the natural and realistic framework 
where to study the stability of the plasma disk. 

2 Two-Dimensional MHD Model for 
an Accretion Disk 

We will now set up the basic 2D MHD configuration equa- 
tions for an axially-symmetric disk around a compact (few 
Solar masses) and strongly magnetized (a dipole-like field 
of about 10^^ gauss) object, i.e. a typical pulsar source. 
In such a system, we neglect the self-gravity of the disk, 
so that the gravitational potential X) in cylindrical coordi- 
nates (r, z, (/)), has the form 



(1) 



where is the mass of the central body, and G the gravi- 
tational constant. It is worth noting that the axial symme- 
try prevents any dependence on the azimuthal angle (j) of 
all the quantities involved in the problem. The continuity 
equation associated to the mass conservation, i.e. 



V • {ev) = , 
admits the following explicit solution 



ev — -VO A -I- eurSd 



(2) 



(3) 



Here e denotes the mass density, v the velocity 
field of the disk, and 0(r ,z) is a function which 
is odd with respect to the z variable, in order 
to deal with a non-zero accretion rate A/d. In 
fact, if H is the half-depth of the disk, Md results 
to be given by [Bisnovatyi-Kogan and Lovelace(2001)[ 



Montani and Benini(2011)| 



Ma = -2TTr [ evrdz = AnQir ,H) > . (4) 

The dynamics of such a stationary, viscoresistive sys- 
tem is then described by the following equation 
Bisnovatyi-Kogan and Lovelace(2001)| 



-f eVx + Vfe 



Stt 



B- 



(5) 



D being the viscosity and I the unit tensor; the magnetic 
field B, because of the axial symmetry, can be expressed 
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through the magnetic flux function = ip (r , z'^'j in the 
general form 

B =-^tP Ae^ + -e^, (6) 
r J. f 

where / = I{ip , z) is a generic function of ip characterizing 
the toroidal component of the field. The similarity of the 
magnetic field and matter flux structure, is due to their 
common divergence-less nature. A further assumption in 
our derivation is that the angular frequency of the disk uj 
has to be a function of "0, i-e. 



UJ = io (ip) . 



(7) 



This result has been firstly demonstrated by Ferraro in 
[Ferraro(1937)| under the hypotheses of an ideal, non re- 
sistive, axisymmetric rotating plasma. It can be shown (see 
the detailed proof discussed in the last section) that this 
result holds even in the presence of small resistivity and 
viscosity. 

2.1 Local perturbative analysis 

In the present analysis, we are interested in the behavior 
of the configuration around a fixed, fiducial value tq of 
the radial variable r. In order to investigate the effects 
induced on the disk profile by the electromagnetic reaction 
of the plasma, we split the mass density and the pressure 
contributions in the zeroth (barred) and the first (hatted) 
orders components, i.e. 



e{ro,z'^) + e, p ^ p{ro, z"^) + p , 



(8) 



respectively. The same way, we express the magnetic sur- 
face function as the sum of a term ipQ, describing the field 
of the star, plus a correction ^pi 



^ = -ipniro) + V'i(ro ,r - tq ,z^) . 



V'l < "00 ■ 



(9) 



The quantities e, p and tpi describe the change induced by 
the currents which rise within the disk. In general, these 
corrections are small in amplitude but with a very short 
scale of variation. Having this scenario in mind, we ad- 
dress the so-called "drift ordering" for the behavior of the 
gradient amplitude, i.e. the first order gradients of the per- 
turbations are of zeroth-order, while the second order ones 
dominate. Accordingly, the profile of the toroidal currents 
rising in the disk, has the form 



As a consequence of the corotation theorem, in the 
present split scheme, we can decompose uj in the sum of the 
Keplerian term ojk and a constant lo'q, respectively given 

by 

w = + tjj'ntpi (12) 



WK = w (V'o) = 



IGM, 



duj 



(13) 



4'o 



This form for uj holds locally, as far as (r — tq) remains a 
sufficiently small quantity, so that the dominant deviation 
from the Keplerian contribution is due to ipi only. 

2.2 Vertical dynamics 



As developed in [Coppi(2005)| 

Coppi and Rousseau(2006)"] , we can distinguish the 
fluid components from the electromagnetic back-reaction. 
The vertical component of eq.([5]) implies that, at the 
zeroth order, the background profile is determined by 
the pure thermostatic equilibrium standing in the disk 
when the vertical gravity (i.e. the Keplerian rotation) is 
sufficiently high to provide a confined thin configuration 



D(z' 



eo(fo) 
e(^o, 0), 



cxp 



z 

IP 



(14) 



Here we have assumed that the background obeys an 
isothermal relation, and that the temperature T may be 
then represented as 



P 



^l±P_^2Kji{f + f), (15) 



Kb being the Boltzmann constant and m the ion mass. On 
the other hand, at the first order, the vertical equilibrium 
is described by the following equation 



dzP + uj'^ze 



1 



A^A^i = . 



(16) 



2.3 Radial equilibrium 

The radial component of eq.([5]) fixes the equilibrium fea- 
tures of the rotating layers of the disk, and can be decom- 
posed into the dominant character of the Keplerian angu- 
lar velocity, plus an equation describing the behavior of the 
deviation, i.e. 



J0 = — VAB 



47r 



j=^(a;*, + «,'„)-i^Aft,(io) 



On the other hand, the azimuthal component of the 
Lorentz force is related to the existence of the function 
I{ip, z) and results to be equal to 



(11) 



2wKro(e -I- ejuj'^^ipi 



— dr 



P + 



{drlplY 



dr^idlipi 



(17) 



Here we neglected the presence of the poloidal currents, 
associated with the azimuthal component of the magnetic 
field, as soon as their contribution resiflts to be of higher 
order. 
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2.4 Azimuthal equation 



By inserting eq.(12| into the ^-component of eq.([5|, we 
easily obtain the relation governing the azimuthal equilib- 
rium; more precisely we get 

OJr. 

1 1 (18) 

= —dr {Dr'cv'odr^) + —dz (Bruj'odzij) + 4 . 
r lUq oJq 

By virtue of the magnetic field form (|6]), of the Lorentz 
force (111, the azimuthal equation stands at tq as 



LOK 



where Dq = lD(?'o) and we neglected the 2-derivatives of the 
viscosity coefficient in view of the hierarchy of the equa- 
tions Montani and Benini(201l")] . This relation accounts 
for the angular momentum transport across the disk, by 
relying on the presence of a viscous feature of the differ- 
ential rotation. Furthermore, this equation establishes a 
link between the radial and vertical velocity fields with the 
poloidal currents present in the configuration. 

A complete analysis of the dynamics is, however, still not 
possible. The system given by Eqs.(16), (171 and (191 is 



not closed, and we have to impose some other condition in 
order to describe the poloidal component of the magnetic 
field / and the function Q. This task will be accomplished 
in the next Section when the balance of the forces acting 
on the electrons will be discussed. 



3 The electron force balance equation 

In the presence of a non-zero resistivity coefficient j], the 
equation accounting for the electron force balance reads as 



E 



AB = T]J. 



(20) 



Since the contribution of the resistive term is expected to 
be relatively small, it turns out as appreciable only in the 
azimuthal component of the equation above. Thus, the 
balance of the Lorentz force has dominant radial and ver- 
tical components, providing the electric field in the form 
predicted by the corotation theorem, i.e. 

E = -- /\ B = -- (driper + dzi^Cz) . (21) 
c c 

Since the axial symmetry requires E^ = 0, the (f)- 
component of the equation above takes the form 



VzBr - Vr-Bz = CrjJ^ , 



(22) 



and, taking into account eq.(lO), it reduces to 

„2 



VrB, 



V,Br = 



rjc 







(23) 



This relation, together with the azimuthal equilibrium 



equation (19), provides a system for the two unknowns 



& and /, when ipi and e are given by the vertical ( 16 1 and 



in eq.(19), we get 



radial (17) equilibria. More precisely, substituting eq.(23) 



2eu 



UJK 



ro Do 



in 



(24) 



+ -. 2~ [^z/ {droi'O + drlpi) - dridzlpl 



(19) 4 Non-linear Configuration 



In order to analyze the dynamics, we need to state some 
relation between the viscosity D and the resistivity rj of the 
plasma. As already done in Montani and Benini(201l")] , 
we assume that the magnetic Prandtl number of the plasma 
is of order one. This in turns implies that we can take the 
following relation 

c"^(e + e) 



Now, inserting eq.(25l into eq.(24), we can get a direct 



relation between the radial velocity Vr, that accounts for 
the accretion on the central object, and the Lorentz force 
F^; more precisely, we get 

Vr = „ i^^^ P^ ■ (26) 
2wK(e + e) 

Substituting Vr by this equation into the electron force bal- 



ance ( 23 ) , we get the following expression for v 
1 



2UJK 



(27) 



In terms of the function Q, the expressions above for Vr 



( 26 1 and Vz ( 27 ) stands as 



dzi>i 



(28) 
.(29) 



The solution of the first equation can be taken as 

e = -^/'^^(^0), (30) 

which, substituted in the second, yields an integro- 
differential relation for F^, i.e. 

dz^pi / dz{drF^) = (c'roV'o + drtpi) F^ - 2wk]DoAi/'i . 



(31) 

Once the behavior of ipi is provided, from the equation 
above we can determine the form of F^{r, z) and eventually 
calculate the (^-component I{'4>, z) of the magnetic field. 
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4.1 Dimensionless variables 



Let us define the dimensionless functions Y, D and P, in 
place of ipi, e and p, i.e. 



Y 



OroVo eo Po 



where 



Po_ 

02 



3e2 



Po = IK^iTtQ , fco = 3^ , v\ 
vi 



2 _ B^o 



(32) 

(33a) 
(33b) 



and B^o ~ Bz{r,z = 0) = drgipo/rQ. We introduced 
the fundamental wavenumber ko that quantifies the typ- 
ical length scale of the perturbation along the radial direc- 
tions. Finally, we rescale the radial variable x = fco(r — ro), 
and the vertical one u = z/Lz assuming the fundamental 
length in the vertical direction to be = JT^H . 



By these definitions, the vertical ( 16 1 and the radial equi- 

(34a) 



libria (17 1 can be restated as 



i? + ^^^r + A,y(i + a,y) + ^a,p = o, (34b) 



where 



A,y = dlY + ^zdlY . (35) 



The two equations above provide a coupled system for P 
and Y once the quantities D and D are assigned; hence 
we are able to determine the disk configuration induced by 



the toroidal currents. Finally, eq.(31) takes the form 

a„r y" duid^A^) = -A,y + (i + a,y) , (36) 

A^ = F^jlijjYj^QkQ being a dimensionless function. 

5 The ring sequence 



Here we will derive a solution to the system (34) 
in the case tz 7^ 0, showing how the formation 
of the crystalline structure for the magnetic field 
and the decomposition of the disk in a ring se- 
quence can happen also in different regimes from those 
addressed in [Coppi(2005) Coppi and Rousseau(2006)"| 
Montani and Benini(2011)| . Let us look for a particular 
solution of the form 

Y (x, V?) = F (m^) sin {ax) , (37a) 

B (x, u^) = F {u^) e^'^'^C'^^) , (37b) 
P (x, 1*2) = M {v?) sin^ {ax) -f N {u^) e''™«('^^) , (37c) 



where a, h are some parameters to be fixed. We assume to 
be very close to the equatorial plane, in order to approxi- 
mate the function 13 to a constant 

|u| < ^ ^ B (u^) = exp {-EzU^) ^ 1 . (38) 



Inserting eqs.(37) in (34), the following system can be eas 
ily obtained 

F" - B^F = , M ^F^ , 



ab ' e. 



(39a) 
(39b) 



where ()' = d/dit, () = d/dv?, and B a positive constant 
x2 - 1 



52 



> 



> 1. 



The solution to system ( 39 ) reads as 
F = Ae-^^, 



N = Cexp 



abC 
i = „ „ ■ exp 



3A(a2 - 1) ' 
abil- bVu^] 



1 - BVu 



6£?A 



3^ (a2 - 1) 



(40) 

(41a) 
(41b) 
(41c) 

, (41d) 



where A and C are integration constants. The finiteness 
of the solutions, leads us to require 

^>0. (42) 
Finally, the total mass density assumes the following form 



B{x,u'^) ^D+^B = 



= 1 



abC 
~2A 



exp 



abe 



3A (a2 - 1) 



U2 - 1 



X exp ( B\fvP- + b cos {ax) 



(43) 



which results to be positive definite if the following condi- 
tion holds: 

The formation of the ring sequence can be seen in Figs.([T]) 
and ([2| representing the density ( 43 ) for a particular choice 
of the model parameters. 



abC 




2A 


exp ^ 



It is worth noting that solutions (41) satisfy also eq.(36) 
as soon as we fix the form of the Lorentz force as the fol- 
lowing 

A^ (x, u^) ^ -F {u^) sin {ax) . (45) 



Finally, from eq.(26), we obtain a non-zero radial velocity 
ADoko e-^^ sin {ax) 



eo 



D 



(46) 



5 



D{x,u ) 0.5 




Figure 1: 3D-plot for the total density D (43|), with 



parameters chosen as A = 1 , C = 
-Vl5,& = -|,£z = 0.05,5 = ^/40. The den- 
sity profile exhibits an modulation along the ra- 
dial direction, while in the vertical direction two 
minima appear. 



1.0 
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1.0 1 
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Figure 2: Level lines for the total density Z)(x,u^) (43), 
with parameters chosen as in fig.(IT|). We note 
that the darkest zones are the minima of density, 
and the ring sequence appears in the equatorial 
plane. 



6 The corotation theorem for a 

weakly resistive and viscous plasma 

In the following, we will show how the corotation theorem 
by Ferraro |Ferraro(1937)| holds under the hypothesis of 
stationarity up to the second order in small dissipative ef- 
fects, poloidal velocity and toroidal component of the mag- 
netic field. 

Let us consider the Maxwell equations 



E+ -vAB, 
c 

VAB= —J, 
c 



and combine them by eliminating J 

^2 



cE ■ 



47r 



r]\/AB-vAB . 



(47) 
(48) 

(49) 



In the case of a steady-state configuration, the vector E 
is irrotational. So, as soon as the poloidal velocity and 
the toroidal magnetic component are small enough, and 
by using eq.([3]) and eq.([6|, eq.(49l reduces to 



VA 



^-Lodr-ip + r]-dzl 
r 



\ 



47r , 



V2.Br) + T] 



47r 



^— + 
r 

l]-drl 

r 



drip 

r 



(50) 

The radial and vertical components of eq.(50), which are of 



first order in the small quantities, are identically satisfied 



as soon as one imposes the electron force balance (22). On 



the other hand, the azimuthal equation can be rewritten 
for convenience in the vector form 



47r / - - \ ^ VI 
(^Vcj A V-0j = VA I ^ — A (/) 



(51) 



The corotation theorem, stating that Vw A Vi/; — 0, can 
be recovered up to the first order because the right hand 
side of the last equation is of second order in 7y, /. Indeed, 
we have to require not only t] and I to be small, but also 
that their gradients must retain the same behavior. It is 
possible to check that this is the case allowed in the solution 



derived above (see the crossmath of eq.(ll) and eq.(45)) 



7 Concluding Remarks 

We analyze an axis-symmetric thin disk configuration 
whose structure is summarized by a mainly rotating plasma 
embedded in a magnetic field endowed with a small toroidal 
component. The consistence of the configuration scheme 
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requires that correspondingly small poloidal matter fluxes 
and dissipative effects must be included in the problem. 
This analysis, which generalizes the original approach by 
Coppi and collaborators, has the merit of providing a co- 
herent description of the vertical and radial dependence of 
the plasma profile. Matter fluxes result to be necessary 
when speaking of a real accreting disk, while the toroidal 
component of the magnetic field is expected to be non zero 
when dealing with the magnetic configuration of an astro- 
physical source. The profile we obtained corresponds to 
have a magnetic Prandtl number of order unity (which is 
a typical value required to phenomenologically address the 
observations |Spruit(2010)| ), and it offers the natural sce- 
nario where a perturbation scheme of the ring sequence 
profile can be implemented. 

This work was developed within the framework of the 
CGW Collaboration (http:// www.cgwcollaboration.it). 
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